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Abstract 



We demonstrate that the relic neutrino asymmetry evolution equation derived 



*^ from the quantum kinetic equations (QKEs) reduces to the Boltzmann limit 

i ' that is dependent only on the instantaneous neutrino distribution functions, 

in the adiabatic limit in conjunction with sufficient damping. An original 
. physical and/or geometrical interpretation of the adiabatic approximation is 

given, which serves as a convenient visual aid to understanding the sharply 
contrasting resonance behaviours exhibited by the neutrino ensemble in op- 
posing collision regimes. We also present a classical analogue for the evolu- 
tion of the difference in the v a and u s distribution functions which, in the 
Boltzmann limit, is akin to the behaviour of the generic reaction A # B with 
equal forward and reverse reaction rate constants. A new characteristic quan- 
tity, the matter and collision-affected mixing angle of the neutrino ensemble, 
is identified here for the first time. The role of collisions is revealed to be 
twofold: (i) to wipe out the inherent oscillations, and (ii) to equilibrate the 
v a and u s distribution functions in the long run. Studies on non-adiabatic 
evolution and its possible relation to rapid oscillations in lepton number gen- 
eration also feature, with the introduction of an adiabaticity parameter for 
collision-affected oscillations. 
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I. INTRODUCTION 



The effects of decohering collisions and coherent flavour oscillations on a neutrino en- 
semble are collectively quantified by the quantum kinetic equations (QKEs) which, in 
recent years, have been frequently applied in the study of neutrino asymmetry (difference 
in neutrino and antineutrino number densities) evolution in the early universe, with appre- 
ciable success IHH- In essence, the simplest scenario involves independent oscillations of an 
active neutrino v a and its antineutrino V a , commencing with equal number densities, to a 
corresponding, initially absent light sterile species, v s and V s , in an environment that alters 
the two sets of oscillation parameters in dissimilar ways. Both v a <-> v s and z/ a <-> V s systems 
evolve subject simultaneously to a biased collision scheme which, in crude terms, is blind to 
the sterile variety. For the appropriate vacuum oscillation parameters, numerical solutions 
to the exact QKEs demonstrate the combined effect to be one that sees the relic a-neutrino 
asymmetry grow to orders of magnitude above the baryon-antibaryon asymmetry |3|-|)|. For 
other relevant works see Ref. 

As much as one would like to use the QKEs directly in all applications, obtaining exact 
numerical solutions remains a computationally intensive task, given the necessity to track 
neutrinos at all momenta where nonzero feedback couples their development. Several ap- 
proximate treatments have been employed in the past which, besides lifting the burden on 
the computer considerably, also on occasion offer valuable analytical insights on the nature 
of the asymmetry evolution PJ^JTI1[]. Two such schemes are the well-established static ap- 
proximation mj, and the closely related adiabatic limit approximation introduced in Ref. 

Beginning with the QKEs, the adiabatic limit approach comprises a set of systematic ap- 
proximations, leading ultimately to an approximate expression for the neutrino asymmetry 
evolution. The extraction procedure, however, was hitherto largely motivated by mathe- 
matical convenience; physical, or at least geometrical, interpretations of the intermediate 
steps and quantities that arose therein were lacking. It is thus our intention in this paper to 
examine the derivation of the adiabatic approximation more closely, and to assign definite 
physical meanings to as many mathematical manoeuvres as possible. These new interpre- 
tations are most useful for visually tracking the evolution of the neutrino ensemble across 
a Mikheyev-Smirnov-Wolfenstein (MSW) resonance [TIJ in different collision regimes, and 



in understanding the nature of non-adiabatic evolution. We shall also demonstrate again, 
this time more assertively, that the adiabatic procedure leads to the elimination of phase 
dependence in the regime of interest: the approximate evolution equations for the neutrino 
distribution functions depend only on the distribution functions themselves, and not on the 
coherence history of the ensemble. In other words, the QKEs yield classical-like Boltzmann 
equations in the adiabatic limit. 

Our second goal in the present work is to draw attention to the similarity between 
some aspects of the QKEs and the behaviours of the more familiar classical linear harmonic 



x The magnitude of the final asymmetry found in Refs. HH was recently questioned in Ref. 
Reference |I will explain the nature of the error in Ref. H . 
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oscillator. In particular, the evolution of the difference between the v a and v s distribution 
functions in the adiabatic limit may be modelled by a damped harmonic oscillator with a 
decaying oscillation midpoint. This analogy turns out to be a very illuminating one that 
illustrates clearly the reason behind the remarkable accuracy and the conceptual correctness 
of the original heuristic static approximation. A new quantity, the neutrino ensemble's 
matter and collision-affected mixing angle, is also established in the course. 

Before proceeding, let us state plainly what we hope to achieve ultimately from these 
abstract analyses: Acquisition of a clear understanding of the physical processes that con- 
stitute the various computationally convenient approximations will leave us better equipped 
to improve on them. To this end, classical analogies are especially useful as visual aids. On 
a grander scale, studies of the QKEs per se may stand to benefit areas beyond relic neu- 
trino asymmetry evolution, most notably transitions in multi-level atomic systems which 
are described by similar equations [0 . 



The structure of this paper is as follows: The exact QKEs are presented in Sec. [TT] for 
the purpose of introducing the nomenclature. Section |TTT| is devoted to the discussion of 
the formal adiabatic procedure, in which we shall also present results from numerically 
integrating the pertinent approximate neutrino asymmetry evolution equation for the first 
time. The classical harmonic oscillator analogy is to be treated in Sec. [TV], while Sec. [V] 
deals with the issue of non-adiabaticity and its possible relation to rapid oscillations in the 
asymmetry evolution. We conclude in Sec. [VI|. 



II. NEUTRINO ASYMMETRY EVOLUTION AND THE QUANTUM KINETIC 

EQUATIONS: NOMENCLATURE 

We shall consider a two-flavour system consisting of an active species v a (where a = e, 
/i or r), and a sterile species v s , where their respective abundances and the ensemble's 
coherence status at momentum p are encoded in the density matrices 

p(p) = ~[Po(p) + p(p)-4 (i) 

The vector P(p) = P x (p)x + P y {p)y + P z (p)z may be interpreted as the "polarisation" and 
a = a x x + <j y y + a z z are the Pauli matrices. In this notation, the v a and v s distribution 
functions at p are respectively 

N a (p) = ^[P (p) + P z (p)]N c ^p,0), 

N S (P) = \[Po(p) - P,(p)]N*(p, 0), (2) 

for which we have chosen the reference distribution function N eq (p, 0) to be of Fermi-Dirac 
(equilibrium) form, 
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with chemical potential [l set to zero at temperature T. The four variables Pq(p) and P(p) 
advance in time according to the quantum kinetic equations (QKEs) 



<9P 

~dt 
dPo 
dt 



V(p) x P(p) - L>(p)[P x (p)x + P y (p)y] 
'iV eq (j»,/i) 1, 



dt ' 



r(p) 



iV ec i(p, o) 



rFo(p) + ^ 



(4) 



where the quantities -D(p) = T{p)/2 and V(p) individually characterise the collision-induced 
decohering and the matter-affected coherent aspects of the ensemble's evolution respec- 
tively.0 Their approximate forms shall be detailed shortly. Note that the equation is 
not exact because the right hand side assumes thermal equilibrium for all species in the 



background plasma, while the v a distribution is taken to be approximately thermal [10 



The properties of the antineutrino ensemble may be parameterised in a similar manner and 
subject to the same QKEs with the appropriate decoherence function and matter potential. 
Henceforth, all quantities pertaining to the V a «-> V s system shall bear an overhead bar. 
The vector V(p) reads |0 



V(p) = (3(p)x + X(p)z, 



(5) 



with 



Kp) 



Am 2 
2p 

Am 2 
2p 



sin 29 , 



[Kp) ~ a (p) ~ cos 2^o 



(6) 



in which Am 2 is the mass-squared difference between the neutrino states, do is the vacuum 
mixing angle, and 



a{p) 

Kp) 



4:((3)^/2G F L^T 3 p 



7T 



''Am 2 



4((3)V2G F A a Ty 



7T- 



'■Am 2 



rn 



(7) 



given that G F is the Fermi constant, m\y the W-boson mass, ( the Riemann zeta function 
and A e ~ 17, A^ T ~ 4.9. The quantity 



L {a) = L a + L e + + L T + 7] = 2L a + L, 
combines the various asymmetries individually defined as 

L, 



(8) 



(9) 



2 It shall be demonstrated later that the system's damping and oscillatory features are in fact 
mutually dependent. 
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where n a = J N a (p)dp, n 7 is the photon number density, and rj is a small term related to 
the cosmological baryon-antibaryon asymmetry The authors of Ref. have called the 
effective total lepton number (for the a-neutrino species), a name we shall also adopt.Q The 



condition A = is identified with a Mikheyev-Smirnov-Wolfenstein (MSW) resonance [11 

The function D(p) is equivalent to half of the total collision rate for u a with momentum 
p, that is 0,0, 

Dip) = \m = l^ly^lT'i 1 - Z M + °( L l)l ( 10 ) 

where (p)o ~ 3.15T is the average momentum for a relativistic Fermi-Dirac distribution 
with zero chemical potential, y e ~ 4, y^ T ~ 2.9 and z e ~ 0.1, z^ T ~ 0.04. 

The corresponding functions V(p) and D(p) for the antineutrino system are obtained 
from their ordinary counterparts by replacing and L a with —L^ and — L a in Eqs. (^) 
and ( |T0| ) respectively. 

We conclude this section by noting that a direct time evolution equation for the neutrino 
asymmetry may be derived from the QKEs together with a + s lepton number conservation 
[1Q[1 , which reads 



dL r 



dt 2n 7 



J (3[P y {p) - P y (p)]N^(p,0)dp. (11) 



As well as serving as the backbone on which to develop useful approximations, this expression 
is useful when numerically solving the QKEs. Although it is redundant, Eq. (|11]) has the 
virtue of tracking the crucial quantity L a without the need to take the difference of two 
large numbers. 



III. ADIABATIC LIMIT 
A. The Boltzmann limit 



The adiabatic limit approximation introduced in of Ref. JTO|0 consists of firstly setting 
the repopulation function in Eq. (|J) to zero, i.e., ^ ~ 0, such that the QKEs simplify to 
the homogeneous equations 



/CP, (12) 





3 The terms "lepton number" and "neutrino asymmetry" are used interchangeably throughout this 
work. 

4 Note that the derivation we now review is very similar to a procedure described in Ref. J14|| , 
though the context is different. 
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where the dependent variables and coefficients are understood to be functions of both time 
and moment um.0 

We solve Eq. ( |T2| ) by establishing a parameter-dependent "instantaneous diagonal basis" 
S g = S q (D, A, (3) = ( qi, q2, <?3 ), onto which we map the vector P from its original "fixed" 
coordinate system S — ( x, y, z) via 



= Q=UP. (13) 

The transformation matrix U and its inverse U~ x diagonalise the matrix K, in Eq. (|i~2|) by 
definition, 

K d = diag(A; 1 , fa, fa) = UKU~\ (14) 

where the eigenvalues fa, fa and fa are roots of the cubic characteristic equation 

k 3 + 2Dk 2 + {D 2 + X 2 + (3 2 )k + (3 2 D = 0, (15) 

the discriminant of which is identically 

A = 4/3 6 - I3 4 D 2 + 12/? 4 A 2 - 20(3 2 D 2 X 2 + 4D 4 A 2 + 12/5 2 A 4 + 8£> 2 A 4 + 4A 6 . (16) 

Under most circumstances (for example, the individual cases of |A| ^> D and \(3\ ^> D), 
the inequality A > holds such that two of the three eigenvalues occur predominantly as a 
complex conjugate pair, which may be conveniently parameterised as 

fa : 2 = —d ± iuj, (17) 

where d and u, both real and positive, are readily interpreted as the effective damping 
factor and effective oscillation frequency respectively. These phenomenological parameters 
arise since the QKEs couple the damped and oscillatory aspects of the time evolution, and 
are to be compared with what could be called the bare damping factor D and the bare 
matter- affected oscillation frequency Q = \J (3 2 + A 2 . The remaining root fa of Eq. (|15|) is 
real and negative, and bears several simple but elucidating relationships to fa and fa, 

h - (18) 
d = D + ^, (19) 

^ = X 2 + (3 2 + faD + ^kl (20) 



5 The approximation ~ requires careful justification, although past numerical evidence has 
always strongly suggested this idealisation to be very reasonable. One can derive in the adiabatic 
limit of the full-fledged QKEs, including the finite repopulation term, results identical to those 
obtained in this section to leading order. See the companion paper Ref. [15|. 
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to be further discussed in Section |TV|. These equations show that k 3 quantifies the "mis- 
alignment" between the effective and bare damping factor, and between the effective and 
bare oscillation frequency. 

The only (so far) unambiguously identified instance for which the discriminant is zero 
or negative (A < 0) occurs when the conditions D 3> |/3| and A 2 ~ ^jj2 are simultaneously 
satisfied, in which case all eigenvalues are real and generally distinct with the exception of 
A = 0, where at least two roots are equal.Q The interested reader is referred to Appendix 
[A| for the mathematical details. For clarity, we shall deal exclusively with off-resonance 
evolution in this subsection, which translates loosely into requiring |A| > D, \/3\, so as to 
ensure the existence of complex conjugate eigenvalues. The study of resonance behaviour 
is deferred to Sec. I111B. 



i.e., where A ~ 

For A 0, the matrix U 



consists of the normalised eigenvectors 



1 

M 



( 1 

A 

a D+kj 



(21) 



in the columns, while the row vectors 



j^i k% kj 



— n 

(3D 11 



constitute the inverse matrix U, with the ith normalisation factor given byQ 

Are* 



(22) 



(23) 



or explicitly, with the aid of Eqs. fll8|) to (|20| ) 

1 



Mi =M 2 



|A| 



2[\ 2 + P 2 + k 3 (D + k 3 



Ma 



\Xka 



D(D + k 3 



(24) 



In this new basis S q , the unit vector q 3 represents the axis about which P precesses. This 
precession axis coincides, by definition, with the matter potential vector V in the absence of 
collisions (i.e., D — 0), residing entirely on the xz-plane when viewed in the fixed coordinate 
system S. Otherwise the alignment is inexact in the general case, where a nonzero D 
generically endows q 3 with a small ^-component evident in the transformation matrix U^ 1 
in Eq. (j2J), that is, 



6 Equations (|T^) to (^0|) are equally valid for cases where A < if we extend the definition of u 
to include imaginary values, although w's physical meaning is then, of course, different. 

7 A minor issue pertains to the definition of the normalisation factor — to be clarified later where 
relevant. 
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q 3 = U x3 1 x + U y3 1 y + U z3 1 z 



l Z3 



(25) 



where 



U 



x'A 



ff 2 



D(D + k 3 



(D + k 3 y + \ 2 + (3 2 

P 2 D(D + fcg) 
\{D + A: 3 ) 2 + A 2 + {3 2 } 2 
(D + h) 2 + A 2 



(L> + k 3 ) 2 + X 2 + (3 2 



D=0 



/3 2 



A 2 + /3 2 



D=0 



D=0 



0. 



A 2 



(-D + k 3 f + A 2 + /3 2 A 2 + /3 2 



(26) 



The .D = (i.e., (73 || V) limit of Eq. (|26"D is to be compared with Eq. (|5]). The variables Q\ 
and Q2 quantify the actual precession in conjunction with their associating complex vectors, 
qi and ^2, that sweep the plane perpendicular to q^^\ Figure [I] shows schematically the 
relationship between the S and S q bases. 

Rewriting Eq. (|1^) in the new instantaneous S q basis, we obtain 



— = /QQ-W— Q, 



(27) 



where the 3x3 matrix U^- exemplifies the system's explicit dependence on ^ and 
and in general contains nonzero off-diagonal entries. Under certain circumstances, these 
time derivatives are inconsequentially small relative to terms in the diagonal matrix /Q (see 
Ref. [H| for the relevant constraints). When these bounds hold, we are entitled to take the 
adiabatic limit, defined by U^j- ~ 0, or equivalently, 



U(D + 5D, A + 5X, (3 + 5/3)U-\D, A, (3) ~ diag(l, 1, 1) 



(28) 



such that the basis S q (D, A, (3) maps directly onto the neighbouring S q (D + 5D, A + <5A, (3 + 
8(3) as the parameters D, X and (3 slowly vary with time. 

The formal solution to the now completely decoupled system of differential equations for 
Q is simply 



Q(t)^diag(e/o fc ^', e Io k ^ d A Q(0), 

and consequently, 

P(t) -ZT^diag ( e Io k ^ dt ', efoWW, X k ^' )dt ') W(0)P(0). 



(29) 



(30) 



Note that this solution is "formal" because the eigenvalues ki depend on the asymmetry L a , 
and hence on P z , through the function A. 



This description of the S q coordinate system is actually technically inaccurate since the instan- 
taneous diagonal basis is not orthogonal, and are complex. However, as a visual aid, it is more 
than adequate. 
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At this point, we alert the reader to a mathematical subtlety. Equations (|29|) and 
therefore (|30"D are in fact poorly defined if the time integration encompasses regions where 
two or more eigenvalues are identical. In this case, the matrix K, has momentarily less than 
three distinct eigenvectors, thereby rendering the transformation matrix U~ l uninvertible. 
Fortunately, such instances, besides their rarity (see Appendix [A]), have virtually no bearings 
on the outcome as we shall see in due course. 

From Eq. fl3"0|), we extract a formal expression for P y (t) in terms of P z (t), 

p , t) „ u^if)^^^ (31) 

v U^{t)J^ l{t ' )dt 'QM+U^{t)Jl k2{t '^ 



Recall from Eq. fllTD that the eigenvalues k\ and k,2 are prevailingly composites of a real 
damping factor d and an imaginary oscillatory component u. Given an ample A such that the 
effective oscillation frequency reads uj 2 ~ A 2 in Eq. fl20|), the real eigenvalue k$ is guaranteed 
to be small according to Eq. (fl8|) , while the phenomenological damping becomes d ~ D 
by Eq. ([OJ). Since k 3 scales with D, a comparatively large d will always quickly reduce 
the corresponding exponentials in Eq. (|3lD to zero relative to the "decay" time of their k 3 
counterpart, wiping out the accompanying rapid oscillations in the process. The condition 
D ^> \f3 1 is a bonus which contributes to accomplishing the said exponential damping at 
an even faster rate over the time scale of the k% decay. Installing the collision dominance 
approximation^ 



in Eq. fl3~H) , we obtain 



; /o*i,a(f)*'_ >0j (32) 



p -« - W) m = j rm (33) 



which, together with its antineutrino analogue, allows us to express the neutrino asymmetry 
evolution equation in Eq. (|TTD in an approximate form, 



dL a 1 
dt 2n 7 



/ {k 3 [N a (p) - N,(p)] - k 3 [N a (p) - N s (p)]}dp. (34) 

Observe that Eq. (|3~4T ) involves only the v a and v s distribution functions at any given time. 
Interestingly, phase dependence has been eliminated by the adiabatic procedure in conjunc- 
tion with appreciable damping (collision dominance); the coherence status of the neutrino 
ensemble has minimal influence the asymmetry's time evolution. Geometrically, the quan- 



tity P y (t) / P z {t) in Eq. (33) is but the instantaneous ratio of the y- and ^-components of the 



9 The resulting approximate evolution equations actually involve the joint action of collisions and 
oscillations. However, we adopt the phrase "collision dominance" because of the crucial role played 
by damping. 
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axis qs in the limit of zero-amplitude precession. We shall henceforth refer to this condition 
as the Boltzmann limit of the QKEs. Note also that Eq. (|34f) has the form of a classical- 
style rate equation [see Eq. ( fj6|) below] with the neutrino and antineutrino transition rates 
computed to be \k%\/2 and \k^\/2 respectively. 

For computational purposes, an auxiliary expression describing sterile neutrino produc- 
tion for each momentum state p may be obtained similarly in the Boltzmann limit JT(| . 



d 
dt 



Ns(p) 
N c *{p, 0) 



N a (p)-N s (p) 



(35) 



to be employed in tracking the quantity N s (p) in Eq. 
the same equation is approximated to be 

N a (p)~N e *(p,ti), 



N e *(p } 0) 

!4T). The v n distribution function in 



(36) 



with the interpretation that the neutrino momentum state p is instantaneously repopulated. 
Thus Eqs. (134]) to (^6)) form a fully serviceable set to be simultaneously solved to give L a as 
a function of time. 



B. Resonance behaviour: the A ~ case 



Reference [H| did not provide a full analysis of collision-affected adiabatic evolution 
through an MSW resonance. Since resonance behaviour is a very important issue in the 
study of lepton asymmetry growth, we now provide a careful treatment of this topic. 

As an incentive to study the vicinity of an MSW resonance, let us observe the solutions 
k to the characteristic equation, Eq. (|1^), at exactly A = 0, 



A=0 



_ D D ± y/D»-W 



(37) 



Evidently, whether the square root evaluates to an imaginary or a real number is conditional 
to the relative sizes of the arguments D and /3. If the former ensues, we are forced to identify 
the root —D with the real eigenvalue k 3 which is now twice the magnitude of the effective 
damping factor d = D/2 such that the off- resonance "damping versus decay" rationale of 
Sec. |111 A| no longer applies. The latter scenario is, contrary to off-resonance behaviour, 
non-oscillatory, and requires a new interpretation. 



1. Case 1: \fi\ ~ D 



There are two distinct situations covered under this heading. The first is where adiabatic 
and collision dominated evolution occurs on either side of the resonance, with the condition 
\/3\ ~ D maintained during the crossing. This amounts to requiring that 



• 2 nfl > o r w i n-24 2 | Am 2 |2/ 

sin 20 ~ 3.6 x 10 y n —. . 



10 



-7 



(38) 
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where y = p/T is a dimensionless quantity not to be confused with y a , and we have used 
Eqs. (H) and fllPP evaluated at resonance, cos29 — a(p) ~ 0. The last approximate inequality 
in Eq. fl38| ) is established by setting y ~ 3.15 (i.e., average momentum), | Am 2 | ~ 1 eV 2 , and 
L( a ) 10 -5 is a typical value for the asymmetry immediately after exponential growth. 

The second scenario is where D is genuinely small enough so that the collision dominance 
approximation cannot be made even outside of the resonance region. This situation typically 
obtains, independently of the vacuum oscillation parameters, at lower temperatures since D 
decreases with temperature as T 5 according to Eq. (pi]). 

It is easy to convince oneself by inspecting the discriminant in Eq. ( |TED that a dominating 
/3 ensures the existence of complex eigenvalues, and thereby preserves the precessive nature 
of the evolution of the polarisation vector P for all A. At A = 0, the solutions to the 
characteristic equation, Eq. ([15]), are 



fcl > 2 L=0,|/J|££> 

fc 3 1 A=0i 



D , V4/^-D 2 

— ± i . 

2 2 

-D, 



(39) 



where the sole real root is always identified as k 3 , which, as mentioned earlier, is clearly 
larger than the real components of k\ and k 2 , i.e., the effective damping factor d = D/2. 
The corresponding transformation matrices evaluate to 



u- 1 



\=0,\B\~D 



7/1 



( o 

i fci 

V2 P 
V ~V2 

fO i/3 
-ip 
VI 





1 k 2 
V2 P 

~72 



4f3 2 -D 2 



4/3 2 - 





D 2 



o , 

°) 

-ik 2 
iki 



_4/32-D 2 
2 

4f3 2 -D 2 





(40) 



/ 



where the precession axis q% now lies in the x-direction, and the vectors q\ and q 2 trace out 
accordingly a surface roughly parallel to the ?/z-plane. 

Let us now track the evolution of P, as depicted in Fig. ^], for the situation where 
collision dominance holds on either side of the resonance crossing. As the matter potential 
vector V, initially at (+/?, 0, +A), is rotated on the xz-plane and decoherence function D 
independently modified, the instantaneous unit vector ^3 moves through a continuum of 
predefined positions in the (+x, —y, +z) block dictated by the transformation matrix 
while P precesses about it. The vector q 3 gradually aligns with the x-axis as we approach 
A = 0, such that the decay of Q 3 becomes increasingly dependent on the parameters that 
govern the evolution of P x . This is manifested in the escalating magnitude of the eigenvalue 
&3, which reaches a zenith of | — D\ when q 3 and x are completely parallel at A = by Eqs. 
([39]) and (|40"|) , and Q% and P x are momentarily equivalent.^ Meanwhile, the two variables of 
interest, P y and P z , become progressively more oscillatory as we cross the resonance region 
for the following reason. 



io A 

more insightful interpretation of k 3 will be provided later in Sec. ffyj. 
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In the adiabatic limit, the precessive behaviour of P is parameterised by the complex 
conjugate variables Qi and Q2, and driven by the imaginary component of the corresponding 
eigenvalues k\ and k 2 . The magnitudes of the oscillatory components in P x ,y, z are therefore 
dependent on two factors: (i) exponential damping of Qi and Q2 through the real part of 
k\ and k 2 , and (ii) the "intrinsic" oscillation amplitudes prescribed by the transformation 
matrix W -1 , or equivalently, the projections of the instantaneous eigenvectors q\ and q\ 
onto the fixed S basis. For the ~ D case, Eq. (RDI) shows that the effective damping 
factor d in k\ and k 2 diminishes from ~ D to an all-time low of D/2. This means that 
the variables Qi and Q2 are no longer preferentially exponentially damped relative to the 
decay of Q3 which now occurs at a comparable rate ~ D. On this account alone, no single 
term in Eq. ( |31]) stands out as the presiding factor in the resonance region. At the same 
time, the evolution of P y and P z is now almost totally described by Q\ and Q2, because of 
a growing alignment between the qiq 2 - and yz-planes. Complete alignment (and therefore 
purely oscillatory behaviour in P y and P z ) is attained at exactly A = as shown in Eq. (pEUp. 

On the other side of the resonance, q 3 resides in the (+x, +y, —z) block while V = 
(+/3, 0, — A) continues to move down the xz-plane (see Fig. 0). The variable Q3 regains 
control as |A| increases, and the system advances in time as discussed in Sec. |III A| . 

The oscillatory terms at resonance play a pivotal role in partially propelling the asym- 
metry growth through the MSW effect (dominant mode for \/3\ 3> D, i.e., negligible collision 
rate). In fact, in regimes where collisions are completely absent (i.e., D — 0), one may 
recover from the approximate solution to the QKEs in Eq. ( |30| ) the common expression for 
the adiabatic MSW survival probability ||10|| . Geometrically, adiabatic \/3\ ~ D evolution 
carries the polarisation vector P from the (+x, —y, +z) to the (+x, +y, —z) block across 
an MSW resonance following the precession axis ^3. But as far as Eq. ([£|) is concerned, 
coherent MSW transitions are evidently dependent on the history of the neutrino ensemble, 
and therefore outside the Boltzmann limit. The approximate neutrino asymmetry evolution 
equation, Eq. (0), thereby collapses in the proximity of a resonance under the condition 
1^1 ~ D. 



2. Case 2: D > \j3\ 

The case of D ^> exhibits vastly contrasting resonance behaviour to the previously 
considered scenario. A genuinely large d ~ D in denominator of k 3 , that is, 

^U^-^^ + W 4 ), (41) 

according to Eqs. ([18]) to (|20|), will ensure the latter's smallness even as A approaches the 
resonance region. 

An immediate consequence is a growing projection of ^3, initially in the (+2, —y, +z) 
block, onto the ?/z-plane, accompanied by an x-component that tends to zero with a decreas- 
ing A in Eq. (p6|). Precession about q 3 is sustained and quantified by the complex functions 
Qi and Q2 until we reach a domain where the imaginary parts of the complex eigenvalues 
totally vanish. We locate this territory by evaluating the discriminant of the characteristic 
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equation in Eq. (|16D with the condition D 3> \(3\ in place (see Appendix [A]). The patch 
turns out to be minuscule, with an upper bound 



A 



4D 2 



(42) 



to be compared with the nominal resonance width |AA| ~ \(3\. Here, the eigenvalues are all 
real and negative, with their relative sizes being the only distinguishing feature. Thus, in 
principle, any one root may be labelled k 3 and the remaining two will automatically satisfy 
Eqs. (|H]) to ( pop with an imaginary u. We restrict the choice by matching the eigenvalues 
at the "complex-real" interface, such that the smallest root is always designated k 3 .f^ The 
others, k\ and k 2 , are dissimilar but concurrently large, thereby providing continued justifi- 
cation for Eq. fl32|) and consequently Eq. ( j33| ) in this regime. As an illustration, we calculate 
the eigenvalues at A = 0, 



k 



iIa=o,.d>| 



k 



2lA=0,B>l 



A: 



3lA=0,B>|/3| 



■D, 





Af3 2 


2 




4/? 2 



—D + 0{(3 2 



D 



+ W 4 )- 



(43) 



Note that the expression for k 3 in fact agrees very well with Eq. (|41|). 

The unit vectors <ji, q 2 and q 3 are similarly matched at the boundary. Transformations 
between the fixed S and instantaneous S q bases are as defined previously in Eqs. fl21"|) and 
([22]), although minimal physical significance may be ascribed to the instantaneous basis 
beyond the textbook description of a 3-space spanned by three real linearly independent 
vectors, where q 3 is, incidentally, most aligned with z. The matrices IA~ X and U evaluated 
at A = are routinely displayed below, 



W 



A=0,D>|, 



u 



A=0,D>|, 



/I 







f\ 







k 2 















§_ 







- 



D 2 -0 2 

(S 2 +k 2 
D 2 -4/3 2 
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(3 2 +k 2 2 
D 2 -4/3 2 

l3 2 +k 2 
D 2 -4f3 2 J 



(44) 



in which the quantities k 2 and k 3 are the relevant eigenvalues listed in Eq. (|43 ), and the ratio 
U~ 3 /U~ 3 l = k 3 /j3 remains unchanged from the A 7^ case. Hence the asymmetry continues 



Once matched, the root k 3 will remain smallest in the "real" region since real eigenvalues must 
never cross therein (except at the boundary) in the limit D S> \/3\, as dictated by a nonzero 
discriminant in Eq. ([l6|). 
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to evolve in the Boltzmann limit at resonance in the D ^> |A|, \(3\ regime, as described 



previously by Eq. (|34|) . 

Observe also in Eq. ([0]) that the vector q 3 is indeed entirely on the yz-plane at A = 0. 
Thus the role of adiabatic evolution in the event of D 3> \(3\ is to bring the polarisation P 
from the (+x, —y, +z) to the (— x, —y, +z) block as we rotate the matter potential vector 
V from (+/3, 0, +A) to (+/3, 0, —A), as shown in Fig. |3|. This is in total contrast with the 
\/3 1 ~ D and, indeed, the D = cases, and entails a new definition for the U~ l and U 
normalisation factor to maintain the mathematical correctness of Eqs. (|29|) and (|3(]), that 
is, 



Mi\o>n = x\k~\ v A2|fci|2 + (/?2 + 1 w|jD + ^ |2 - (45) 

Compared with Eq. (p^), this new definition guarantees that the x-component of q 3 , U~ 3 , 
(and not U~ 3 and U~ 3 , as in the ~ .D case) flips sign across an MSW resonance. 



3. Summary 

In this subsection, the following conclusions have been reached pertaining to evolution 
through an MSW resonance: 

1. Where the condition D ^> \j3\ is satisfied, asymmetry evolution (which is just a special 
case of distribution function evolution) continues to be well described by the Boltzmann 
limit. Sufficient damping is in place to remove the coherence history of the ensemble 
as a necessary and independent dynamical variable. 

2. The case of \/3\ ~ D displays completely different resonance behaviour since the damp- 
ing versus decay rationale supplied by |Re(fci j2 )| 3> |A; 3 | does not hold. The MSW effect 
now dominates the dynamics of the system in a way that cannot be described by a 
Boltzmann limit. 



C. Static approximation 

We present here a brief account of the static approximation MM , partly to demonstrate 



that Eq. fl34|) may be re-derived within this framework, but also to introduce some termi- 
nology for later use. The derivation given here is a development on that in Ref. Mj. 

The static approximation begins with the following observations @||: There are two 
means by which the lepton number of the universe may be modified: (i) oscillations between 
collisions, since matter effects are dissimilar for neutrinos and antineutrinos, and (ii) the col- 
lisions themselves which deplete v a and V a at different rates through physically "measuring" 
the respective v a <-> v s and V a <-> V s adiabatic matter-affected oscillation probabilities. 

This vision of the role of collisions is essentially a wave-function collapse hypothesis (pro- 
jection postulate) that is evidently closely related to the collisional decoherence rigorously 
quantified through the D function in the QKEs. Exploitation of this picture has, in the 
past, led to the acquisition of much physical insight (see Refs. [0,0 )• There is a significant 
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literature on the relationship between the collapse hypothesis and systematic treatments of 
quantal decoherence for open systems The general conclusion seems to be that the 

former often leads to dynamical behaviour close to that derived from the formal "master 
equation" or QKE approach. Using this type of picture for active-sterile neutrino oscilla- 
tions, one can derive a lepton number evolution equation (static approximation) similar to 
that obtained in the adiabatic limit of the QKEs in the collision dominance regime |||K|. 

The approximation we shall discuss here neglects mechanism (i) above. Mode (ii), gen- 
erally predominant in the high temperature regime, is described by the rate equation 



dt 



The reaction rate 



1 



[-T(u a -> u s , p)N a (p) + T(u s -> u a , p)N s (p) 

+ T(u a -> u s , p)N a (p) - F(u s -> u a , p)N s (p) dp. 



r(u a -> u s , p) 



(P(u a -> U s , p)), 



(46) 



(47) 



is the product of half of the collision rate T{p), and the probability that an initially active 
neutrino at time t' would collapse to the sterile eigenstate at the time of collision, t, averaged 
over the neutrino ensemble.^ In the adiabatic limit, the latter is given by 



(P(u n 



S,P)) ~ 2 2\n(t)Q(t'] 



+ n(t) n(t>) 



n(t")df 



~-^(c2H m (t)e2UJt')-\ .s2^,(/-).s2^,f/-')(-c)s 



n(t")dt" 



(48) 



where Q = ^X 2 + j3' 2 is the matter-affected oscillation frequency, c29 m = cos 29 m and s29 m = 
sin 29 m quantify the matter-affected mixing angle, the subscript "ens" denotes a type of 



ensemble average, and all quantities are functions of momentum p. It follows from Eq. ( 48|) 
that V(u a — > u s , p) = T(u s — > u a , p), and likewise for their antineutrino counterparts. 

The ensemble average phase (• • -) ens is computed assuming that collisions disrupt coher- 
ent evolution by resetting the phase term in the brackets to zero. Consider the u a ensemble 
at time t. The fraction that has survived resetting since time t' is postulated to be 



z(t, t') 



exp 



D(t")dt" 



(49) 



Thus the portion of neutrinos that is reset between t' and if + dt' and then propagates freely 
to time t without further encounters is simply 



dz(t, t') = D{t') exp 



D(t")dt" 



dt' , 



(50) 



12 The factor of 1/2 arises from the fact that D equals T/2 rather than T in the QKEs. The 
discovery paper |H| omitted this factor. 
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contributing the phase c29 m (t)c29 m (t') + s29 m {t)s29 m {t') cos $ Q(t")dt" to the ensemble 
at time t. The ensemble average phase is a weighted sum of these contributions originating 
in a time interval extending from t' = ti in the past, to the present t' = t: 



c29 m (t)c29 m (t') + s29 m (t)s29 m (t') cos 
t'=t 



Vt(t")dt" 



D 



t'=ti 



1 ' \n{t)n{t') + n{t)Q{t') 



n(t")dt" 



dt' 



m [t c _ f t D(ndt ,, D(t>)\(t>) df , | ip{t) 



Q(t) 



2n(t) 



n(t>) dt + cx ; 



(51) 



where cc. denotes complex conjugate. Integrating by parts, the first integral evaluates to 



^ t D(t')e-ft D ^ dt "^- 



m 4 



■dt' 



n(t) n(u) n(t) Ju 



\\t) - P t D(t")dt" \{t) X(U) X(t) [* f D{t n )difl 

mt)~ e ! 1 ' 

n 2 (t)' 



dt' 



dt' 



(52) 



The last approximate equality arises from the static approximation which assumes negligible 
dependence on the time derivatives of D and Q. The integration limit t — ti — * oo has also 
been utilised to eliminate the definite integral exp[— J* t * D{t")dt"] in Eq. (p^). 
Applying the same procedure on the second term in Eq. (|51|), we obtain 



[\d + in)e-S>wW\ , ^ dt' 

Ju ' ft(t') D(t') + itt(t ') 



n(t) 

P\t) D{t) 



Vt 2 {t)D{t) + iVt{t) 



Jl (D+ift)dt" P(t) (3(tj 



n(t) fife) d(u) + in(u) 



P(t) /"* - f*(D+ai)dtf' d_ 
Q{t) Ju 6 ' dt' 



m D(t') 
tt(t') D(t') + ifi(tf) 



/3 2 (t) D(t) 



n 2 (t) D(t) + in(t) , 

and similarly for its complex conjugate. Thus Eqs. 

(P(u a -> z/ s , p)) ~ -- ^(•••) ens - 

and the original rate equation in Eq. ((46|) becomes 
dL r 



(53) 



and ([53p combine to give 



2D 2 (t)+Q 2 {t)' 



(54) 



dt 



_L | ^static^^) _ ^ a(p)] _ ^aop7 a (p) _ iV s (p)]} rfp, (55) 



with 
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and a similarly defined antineutrino counterpart /^ tatlc evaluated at time t. Equation (|56|) is 
to be compared to the exact eigenvalue k 3 defined in Eqs. (|Tj|) to (p0|), and to its approximate 
form in the limit D, |A| ^> displayed in Eq. (PD- The exact &3 agrees with ^ tatlc when 
the small k 3 contributions in Eqs. fll9|) and (|20|) are neglected. This clarifies a point of 
confusion raised at the beginning of Sec. II C in Ref. (TIJ that was later resolved, for the 
first time, in Ref. [ITS]. 



D. Numerical results 

For the purpose of book-keeping, we include here a report on the results of numerically 
integrating Eqs. (|34D and (|35D together with the approximation of instantaneous repopula- 
tion given by Eq. fl36|). Calculations are performed on several sets of v a <-> za, oscillation 
parameters Am 2 and sin 2 2# for two choices of k 3 : 

1. Exact fc 3 . The real root /c 3 of the characteristic equation, Eq. (|15l), is computed by 
iterations at each time step for all neutrino momentum bins. In the event of three real 
solutions, we invariably choose the one smallest in magnitude as justified in Sec. |111 B 
This calculation, codenamed exactk3, has not been previously attempted. 



2. Static k 3 . In the statick3 code, we adopt the heuristically derived kf 3 " 110 of Eq. (|56|) 
which, in the past, has very successfully generated neutrino asymmetry growths that 
closely mimic solutions of the exact QKEs ||. The present independent calculation 
incorporates a distinct decoherence function D(p) for the antineutrino ensemble that 
was formerly taken to be identical to its neutrino counterpart. 

Results from exactk3 and statick3 for two representative sets of Am 2 and sin 2 29q are 
shown in Figs. |] and |5|. Solutions to the exact QKEs (qke) for the same oscillation pa- 
rameters, computed independently for this work, are also presented in the same figures for 
the purpose of comparison as well as to demonstrate that the QKEs are indeed numerically 
tractable. 

In all three codes, time integration is achieved using the fifth order Cash-Karp Runge- 
Kutta method with an embedded fourth order formula for adaptive stepsize control |17| . The 



neutrino momentum distribution is discretised on a logarithmically spaced mesh following 
Ref. H], and a summation over all momentum bins is performed at each time step. For 
extra numerical stability, the redundant Eq. ([□]) is also built into the qke code, to be 
simultaneously solved with Eq. This feature essentially serves as a safety net that 

insures against errors arising from taking the difference between two large numbers. 

It is clear from Fig. |4] and |5| that for all practical purposes, exactk3 and statick3 
produce indistinguishable results (and we add in parentheses in conclusion that the former 
code is time-consuming to execute). 
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IV. THE CLASSICAL OSCILLATOR 



Much of the discussion in the preceding section, particularly on the meaning of the 
eigenvalues, was conducted in jargons borrowed freely from the damped simple harmonic 
oscillator. We now present a full classical analogue that models the individual time devel- 
opment of the variables P x , P y and P z . 

Consider a classical system exhibiting damped simple harmonic motion about a midpoint 
xq that simultaneously "decays" with time, such as portrayed in Fig. |B| and by the following 
second and first order ordinary differential equations, 

d 2 (x-X ) n d(x — X ) 2 f \ r, ,rrr\ 

df2 + 27 dt + rf(x ~ Xq) = 0, (57) 

and 

w = (58) 

The quantity £ is recognised as a decay constant, uq is the natural frequency of the oscillator, 
and 7 is the damping factor that couples to the oscillator's velocity to effect a dissipative 
force. 

The variable x may be solved for exactly by inflating Eqs. (|57|) and (|58D conjointly to 
produce a third order ordinary differential equation 

^ + (27 + ^ +(27f + ^)| + ^^0. (59) 
Returning briefly to the QKEs (with negligible repopulation function) in Eq. (12), we observe 



that each component of the polarisation vector P may alternatively be separately solved by 
expanding the said system of homogeneous equations into three mutually independent third 
order differential equations, one per variable. The inflated differential equations are 

03p a2p np 

w + 2D ~w + {D2 + x2+p2) ^ + p2DP = °' (60) 

for the case of time-independent parameters D, A and (3. 

Suppose x is the parallel of P z , i.e., the normalised difference between the u a and v s 
distribution functions.0 Comparison of Eqs. (|59|) and (|60|) immediately leads to the following 
identifications, 

1 = D ~% ( 61 ) 
v 2 = (D-O 2 + \ 2 + (3 2 , (62) 

i = — 2 -. (63) 



13 We shall restrict our discussions to the evolution of P z since it is the only variable that represents 
a tangible quantity, although the analogy is equally applicable to P x and P y . 
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With some minor algebraic manipulation, the real quantity £ can be shown to satisfy the 
characteristic equation in Eq. (15). Thus £ is identically — k 3 , as the reader would expect, 
and its role is to equalise the neutrino distribution functions N a and N s , such as in a reaction 
A B where the forward and reverse reaction rate constants are both \k 3 \ /2. From this 
perspective, the Boltzmann limit simply consists of approximating P z (t) as some midpoint 
Pzo(t) [and similarly for P y (t)], which suffices if the oscillatory component P z (t) — P z o(t) is 
negligible relative to P z o(t) through collision- induced damping and/or matter suppression. 

Furthermore, suppose we wish to extend the definition of the matter-affected mixing 
angle 



sin 2 29 r 



ff 2 



\ 2 + (3 2 ' 



(64) 



to include collision effects by recognising that the denominator of the above is simply the 
matter-affected oscillation frequency. Then replacement with the natural frequency v$ im- 
mediately leads to an effective matter and collision- affected mixing parameter 



(3 2 

sin 2 29 m n = 



2 • 



(65) 



which, together with Eq. (|63|), automatically grants k 3 a most intuitive interpretation; 
the reaction constant k 3 reflects on the individual neutrino state's ability to mix (through 
sin 2 29 m> £)) and to collide (through D). This is in complete agreement with heuristic deriva- 
tions. A quick juxtaposition of the exact k 3 obtained from Eqs. (|62|) and (|63D , 



and the heuristic static k 3 , 



-Dsm 2 29 m n 



(3 2 D 



(D - i) 2 + X 2 + (3 2 '' 



r,static 
K 3 



(3 2 D 



D 2 + A 2 + (3 2 ' 



(66) 



(67) 



exemplifies the static approximation's remarkable accuracy in the Boltzmann limit of the 
neutrino asymmetry evolution. Indeed, this definition of the matter and collision-affected 
mixing angle arises naturally from the QKEs with 



sin 2 29 m>D 
cos 2 29 m n 



+ 



D 



h 

D 



(68) 



that are quite visibly related to the projections of the precession axis q 3 in the fixed S basis, 
as shown schematically in Fig. |7|. We refer the reader to Appendix [B] for a more detailed 
discussion. 



The damping term 7 in Eq. (61) is equivalent to d, the real part of the complex eigenvalues 
as defined in Eq. QT7D, and the corollary 



r = 2 7 + £, 



(69) 
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illustrates lucidly the dual role played by collisions: (i) 27 rapidly damps the oscillations, 
and (ii) £ drives the reaction v a v s to equilibrium in the long run. The second order 
differential equation, Eq. (|57|), evokes oscillations at a damping-affected frequency 



222 
v =1/5-7 



X* + p2 + k 3 D + ^kl, (70) 



that is identified with uj of Eq. ([1TD , the imaginary component of complex eigenvalue. Thus 
collisions seem to modify the oscillation frequency in two ways by, (i) contributing to the 
effective mass to produce a natural frequency that is larger than the D-free frequency evident 
in Eq. (|62"D , and (ii) reducing the natural frequency through the damping term 7, a la classical 
oscillators by Eq. ([70|). The end result, however, is that the two effects seem to negate each 
other to some extent. This is as yet not very well understood. 

V. NON-ADIABATICITY 

Significant variations in some or all of the parameters D, A and f3 over a characteristic 
time scale St ~ (D 2 + A 2 + f3 2 )~* generally result in non-unitary mapping between adjacent 
instantaneous diagonal bases S q (D, A, (3) and S q (D + ^B-St, A + ^St, j3 + §f<^), i.e., 

3D 3\ 38 31 J 1 

U{D + —St, A + —St, (3 + -£st)U-\D, A, (3) ~ diag(l, 1, 1) - U— St, (71) 

where 

3D 3\ 33 3 

U{D + —St, A + —St, (3 + -£st) ~ U(D, A, (3) + -U{D, A, (3)St, (72) 

and we have used ^{UU^ 1 ) = ^-U^ 1 1 = 0.0 This corresponds to the non-adiabatic 

regime where the computationally convenient decoupling of the evolution of Qi, Q2 and 
Q3 is invalid amidst sizeable off-diagonal entries in the matrix ti^ that essentially re- 
couple the named variables in Eq. (|27|). The effects of these off-diagonal terms are most 
discernible near the resonance region where the characteristic time St is a maximum. Full 
analytical quantification of their contribution generally requires one to inflate the homoge- 
neous QKEs in Eq. (|12|) into three independent third order ordinary differential equations, 
which, unfortunately, are not readily soluble even for a simple linear A(t) profile, let alone 
three time-dependent parameters D, A and (3. However, the generic role played by rapidly 
changing parameters may be understood through consideration of the following simplified 
situation. 



14 Predictably, the characteristic time 5t ~ — (D 2 + A 2 + j3 2 ) 2 is the local natural oscillation 



length of the neutrino system as defined in Sec. |V], to be further discussed later in Sec. V B 
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A. A toy model 



Suppose the neutrino ensemble evolves adiabatically up to time t c , at which A = A c , 
D = D c and (3 = /3 C , such that 



P a (to) ~ U- 3 l (t , \ c )e^ k ^ dt 'Q 3 (0), 



(73) 



with a = x, y, z, and substantial damping is implicit. Resonance crossing is instigated 
through the instantaneous switching of A c to — A c at time t c , after which the system continues 
to undergo adiabatic evolution to time t, that is, 



Pa(f)^£Z#(*) e 



-i (A JlMt')dt' 



UiP\f>ci ~X c )Pf}{t c 



In the D, |A| ^> limit, the term J2a[' " "] evaluates explicitly to 



y^Mipjtc, ~^c)M/33(tc, A c 

& 

where we have used 



V2/3 C X C V2(3 C X C 



0, 



D* + A 



2 - 



£> 2 + A 



2 ■ 

c 



£> 2 + A 



2 ' 



u- 1 



I _L A l a 

1 V2\X\ V2\X\ 



DM 



§X 
£> 2 +A 2 
BD 



\ 



V2 V2 D 2 +\ 2 

i i(3 D+i\X\ i§_ D-i\X\ i / 

\ V2D 2 +\ 2 V2~D 2 +\ 2 1 I 



+ 0((3 2 



and 



U\ 



n>\&\ 



I J_^_ J_ i/3 D+i|A] \ 

' >/2 |A| V2 V2D 2 +\ 2 \ 

J_2l L »/3 g-ijAJ 

V^JAI V2 V2D 2 +\ 2 

. /3A /3D -| 

V D 2 +A 2 D 2 +A 2 - 1 



so that to the lowest order in (3 

p,{t) - 



_ 2/3 e A c _ D(t')dt'+ f c k 3 (t')dt' 
D? + xf 



sm 



|A(t')K 



Q 3 (0), 



P 2 (t)~e/o fc 3(*'K'g 3(0) + 0(/3 2 ) _ 



(74) 



(75) 



(76) 



(77) 



(78) 



The final step of letting A c approach the vicinity of the resonance, A c — > D c , where | AA| ~ D 
is to be interpreted as the effective resonance width (see Appendix |B|), completes this toy 
model, 



Py(t)\ 



(3D Pc -f* D[t')+k 3 {t')dt' 
D 2 + X 2 D c 6 



sm 



\X(t')\dt' 



P,(t). 



(79) 
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Comparing Eqs. (|33"D and (Jfflj), we find that the first term in the latter is simply the adiabatic 
result in the D ^> limit. The remaining term quantifies oscillations induced by the sudden 
change in A near the resonance (the connection between oscillations and non-adiabatic effects 
was first discussed in Ref. 0). These 0{f3) oscillations, although eventually exponentially 
damped over a time scale of ~ 1/-D, are most visible immediately after resonance crossing, 
t ~ t c , as the two terms in Eq. ( |79| ) may then be comparable in size. 



B. A more rigorous treatment 



The aforementioned post-resonance oscillations may be shown to arise formally from a 

at ' at miU at ■ 



more exact approach that incorporates the actual, finite time derivatives ^r, ^7 and — 



Following Ref. H], we institute the complex variable 



P = P X + iPy, 



(80) 



together with its complex conjugate P*, which, from Eq. QI2"D, advances in time according 
to the equation 



OP 
~dt 



and similarly for P*. 

The formal solution to Eq. 



P{t) = P*(t) 



-- {-D + iX)P - i(3P z , 
is given by 

-S>-^ dt 'P(0) -i t e~ $>~'^ dt " $P z dt' 



(81) 



-i f e~ f>-W 'pp z dt' 



(82) 

where we have used the initial conditions P x (0) ~ P y {0) — to establish the last approximate 
equality. Note that since we are working in the D ^> limit, even if these initial values are 
not exactly vanishing, sufficient exponential damping will nonetheless obliterate the term 
proportional to P{0) in time as discussed earlier in Sec. [Ill A . 
Applying integration by parts on Eq. (|8"2"D, we obtain 



P(t) ~ -i [\d - t\) e - ti {D - iX)dt 



„ (3P Z 



D-iX 



dt' 



i(3P z 



D-iX 



+ e 



> *PP Z 



t'=t 



D-iX 



+ i 



t'=Q 



e J * 



f*,(D-iX)dt" d ( PPz 



dt' \ D — iX 



dt'. (83) 



Substantial exponential damping again wipes out the second term, such that P(t) and P*{t) 
amalgamate to give 

rt p / rl 



P(t) ~P*(t) 

(3D 



D 2 + X 2 



VD 2 + X 2 



Re(S)e-/> d '"sin (J Xdt" 
+ Im(H)e~ ft Ddt " cos 



Xdt" 



dt', (84) 
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where Re(H) and Im(H) denote respectively the real and imaginary parts of 



iD - X 
VD 2 + X 2 



'ldp 
pdt' 



D-iX 



'dD .d\\ 



dt\ 
dt' 



(85) 



Evidently, the first term in Eq. (|S¥) is simply the adiabatic Boltzmann result in the D 3> \(3\ 
limit. The accompanying integral is a close analogue of the 0{(3) oscillatory term found in 



the toy model in Eq. (|79f) , and may be interpreted as a continuous sum of perpetuating 
oscillations induced by sizeable finite time derivatives |£ and |j| at time t! over the 
history of the evolution (recall that the change in A is instantaneous in the toy model). 
Observe also that in the D 3> \(3\ limit, the quantity ^ = f3P y is a negligible 0{f3 2 ) term 



that may be readily verified by iterating Eq. (84|). 

Thus what remains in Eq. ( |55| ) to first order in (3 is in fact equivalent to elements in the 

that is [PHI , 



matrix 



in Eq. 
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dU 
~dt 
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(86) 



(87) 



that were previously discarded in the adiabatic limit. 

The integrand in Eq. (153) is largest at A ~ 0, i.e., in the proximity of a resonance. In the 



context of lepton number generation, we expect the term % to dominate over 4sr and ^ m 



the region of exponential growth (i.e., 
is also correspondingly rapid. Hence, Eq. 



H to dominate over ^ and |j| " 
a ), where, by Eq. (I7|), the rate of change of A 
51) reduces to 
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from which we identify 1/|S| as the physical resonance width, given |AA| ~ D is the collision- 
affected resonance width in phase space. The corresponding characteristic time of the phase 
factor exp[— ${D — iX)dt"] is 
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VD 2 + X 2 
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D'' 



(89) 



to be interpreted in the D ^> \/3\ limit as the local natural oscillation length of the neu- 
trino system. Thus the degree of non-adiabaticity may be quantified by a collision- affected 
adiabaticity parameter 



T 



D 2 



dX 



dt 



(90) 



A~0 



defined as a comparison between the physical resonance width and the natural oscillation 
length of the neutrino system at A ~ 0, such that the condition T C 1 denotes an ideal 
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adiabatic process, and vice versa. The role of T is illustrated in Fig. || which shows P y (t) as 
a function of A(t), where A(t) varies linearly with time from +A at t = [initial conditions: 
P x (0) ~ P y (0) — 0, and P z (0) — 1 ] to —A, for six different matter density gradients 
The parameters D and /3 are kept constant for simplicity, and P z (t) ~ 1 throughout the 
evolution. 

Observe in Fig. § that a sufficiently adiabatic process always returns a negative P y . Sub- 
stantial non-adiabatic oscillations, however, can carry P y periodically across the zero mark 
immediately after resonance crossing. From the perspective of neutrino asymmetry evolu- 
tion, occurrence around the bulk of the momentum distribution could conceivably alternate 
the sign of the rate in Eq. @§ in a cyclic manner, leading possibly to oscillations in 
the integrated variable L a in the region of exponential growth. This effect has been studied 
numerically in Ref. |3J] (see also Ref. ||18|| ). Although a regime of rapid oscillatory lepton 
number evolution has yet to be completely confirmed [|J, non-adiabatic effects do indicate 
that their existence is likely. This reinforces conclusions first reached in Ref. ||. 

VI. CONCLUSION 

The Boltzmann limit of the relic neutrino asymmetry evolution represents a phase in 
which the rate of the said evolution is dependent only on the instantaneous neutrino and 
antineutrino distribution functions. An associated evolution equation may be extracted from 
the exact QKEs in the adiabatic limit where the matter potential and collision rate both vary 
slowly with time, in alliance with adequate collision-induced damping which serves to wipe 
out the coherence history of the system. In the course of the derivation, we have ascribed 
precise physical and/or geometrical meanings for the adiabatic approximation and the many 
quantities arising therein that were previously lacking. These original interpretations allow 
one to easily visualise the sharply contrasting behaviours exhibited by the neutrino ensemble 
across an MSW resonance in the collision-dominated and the coherent oscillatory regimes 
respectively. 

The time development of the individual neutrino and antineutrino ensembles, in partic- 
ular the difference in v a and v s distribution functions, is mimicked by a damped classical 
system oscillating about a midpoint which contemporaneously "decays" with time. We iden- 
tify this decay constant with the "rate constant" that couples to the evolution equation for 
(N a — N s ). Thus the Boltzmann limit is revealed to consist of approximating the intrinsi- 
cally oscillatory quantity (N a — N s ) to some stable average (N a — N s ) which evolves in time 
in the same manner as (Na — Ng) in the generic reaction A # B with identical forward 
and reverse rate constants. The approximation's validity is guaranteed for small-amplitude 
oscillations through substantial damping and/or matter suppression. The pertinent reaction 
rate constant is essentially a product of the decoherence rate and the neutrino ensemble's 
matter and collision-affected mixing angle, the latter of which is a new quantity identified 
here for the first time. Collisions therefore play two roles: (i) to quickly damp the oscilla- 
tions, and (ii) to drive the "static" reaction v a ^ v s to equilibrium over an extended period 
of time. 

Significant time variations in the damping and mixing parameters, i.e., non-adiabatic 
evolution, are shown to induce comparatively large amplitude oscillations in the system im- 
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mediately after resonance crossing. The degree of adiabaticity is quantified by an eponymous 
parameter first introduced in this work, whose role in collision-affected neutrino oscillation 
dynamics parallels that of its more common MSW-style counterpart in completely coherent 
scenarios. Substantial oscillations in the bulk of the momentum distribution may lead to 
periodic alternation in the sign of the quantity and are thus a prime suspect for the 
generation of rapid oscillations observed by others in the course of the asymmetry evolution. 

Ultimately we would like to improve on the adiabatic Boltzmann approximation, and 
results from the present work have put us on better grounds for its accomplishment. The 
classical oscillator analogy, for instance, seems to suggest a procedure by which to correct for 
phase dependence in the asymmetry evolution. This and other avenues await to be explored 
in the future. 
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APPENDIX A: CUBIC POLYNOMIAL 



We present further information on the derivation of Eqs. (|T^) to fl2"0|), and other conditions 
arising from the characteristic equation in Eq. (fL5]) scattered throughout Sec. III . 



The cubic polynomial under investigation is the characteristic equation of the 3x3 
matrix /C of Eq. fll2|) , 

x 3 + 2Dx 2 + (D 2 + A 2 + f3 2 )x + (3 2 D = 0, (Af ) 

where the coefficients a% = 2D, a 2 = D 2 + A 2 + f3 2 , and a 3 = (3 2 D are real and positive. We 
seek three solutions k\, k 2 and k 3 , some of which may be complex, to Eq. (|A1|) , that is, 



(x — k\)(x — k 2 )(x — k 3 ) = 
x 3 - (kx + k 2 + k 3 )x 2 + {k x k 2 + k 2 k 3 + k 3 h)x - k x k 2 k z = 0, (A2) 

with 

k 1 + k 2 + k 3 = -2D, 
hk 2 + k 2 k 3 + k 3 h = D 2 + A 2 + (3 2 , 

k x k 2 k 3 = -I3 2 D. (A3) 

Equation ( |A3| ) shows plainly that any complex roots must occur as a conjugate pair, and 
the remaining one is necessarily real. For this we introduce the following parameterisation, 

fc 1)2 = -d±iu, (A4) 
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where the quantities d and u are defined to be real and positive. Note that the definition of 
u} may be extended to include imaginary numbers, in which case k\ and k 2 are simply two 
distinct, real roots. Simple algebraic manipulation of Eqs. (|A3|) and (|A4|) leads to 



f3 2 D 



d = D + 

.2 



d 2 + OJ 2 



OJ 



2 ' 

A 2 + (3 2 + k 3 D 



-k 



3- 



(A5) 
(A6) 
(A7) 

(A8) 



which are identically Eqs. fll8|) to (|20"D in the main text. The quantity 

A = -(A; 1 -A; 2 ) 2 (A; 2 -A;3) 2 (A;3-A; 1 ) 2 , 
known as the discriminant, thus characterises the nature of the three roots: 

1. A > 0. One root is real and the other two a complex conjugate pair. 

2. A = 0. All roots are real and at least two are equal. 

3. A < 0. All roots are real and distinct. 



These shall be labelled conditions 1, 2 and 3 hereafter. The general form of Eq. ( |A8| ) in 
terms of the coefficients a 1; a 2 and a 3 of the cubic polynomial may be found in any standard 
mathematical handbook, that is, 



with 



Q 



R 



A = 4Q 3 + 27i? 2 



3a 2 — a\ 



9a±a2 — 27a3 — 2a\ 
27 



(A9) 



(A10) 



Thus the discriminant of Eq. (|A1|) is equivalently 

A = 4/5 6 - ^D 2 + 12/? 4 A 2 - 20(3 2 D 2 X 2 + 4£> 4 A 2 + 12/? 2 A 4 + 8£> 2 A 4 + 4A 6 . (All) 

The reader can easily verify that condition 1 eventuates if the requirement |A| ^> D is met 
for all (3, and similarly for the case ^> D with a variable A. The D ^> |A|, \j3\ affair is 
a trifle more intricate. We begin by keeping only terms up to order (f3/D) A and (X/D) 4 in 
Eq. ( [A 11| ). Then condition 1 is attained if the inequality 
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> o, 



(A12) 



holds. Given that the above quadratic has solutions at 
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the inequality in Eq. ( A12 ) may be similarly expressed in terms of two conditions 



^ >( l o 
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(A14) 



Naturally, the second statement is never true. Hence, in the limit D ^> |A|, \/3\, the discrim- 
inant is positive only for 



,2 ? 

A > 



AD 



2 ■ 



(A15) 



as reported in Eq. fl4"2"|) in Sec. UTE. 



APPENDIX B: MATTER AND COLLISION- AFFECTED MIXING ANGLE 



We begin our discussion by putting forward a conjecture, that the mixing angle cos 26* m ,D, 
collision-affected or otherwise, is the normalised ^-component of the precession axis q% (see 
Fig. |l] for orientation), i.e., 



cos 2 29 m n 



;bi) 



where U 1 is the transformation matrix given by Eqs. ([£]]) and fl2"3"|). It follows that the 
remainder 



sin 2 29 m n = 



+ 



(B2) 



is generically projected onto the xy-pl&ne. In the limit D = 0, the unit vector q 3 is parallel 
to the matter potential vector V with no projection on the y-axis. Equations ( [BID and ( |B2"1) 
thereby reduce to the household expressions readily obtainable from Eq. (|5]). 

The general collision-affected case is somewhat more convoluted. Installing the matrix 
elements from Eqs. ( PTj ) and (|2lf), Eq. ( |B2"D evaluates explicitly to 



(B3) 



sin 2 26 m n 



k\ [A 2 + {D + k 3 



\*kl + (F + kl)(D + k 3 y 
Equations QTSJ) to ( pO"D combine to produce 
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which permits us to make the substitution k 3 (D + k 3 ) 2 = 



-[I3 2 D + k 3 (X 2 + /3 2 )] to arrive at 



siu 26 m , D = --= {D + h)2 + X2 + ^ . (B5) 

We identify the denominator as, in the classical harmonic oscillator vernacular, the "natural 
frequency" v% given by Eq. (with k 3 = — ^). Thus Eq. (|E5|), with its lofty QKE origin, 
is entirely equivalent to the heuristically derived matter and collision-affected mixing angle 
in Eq. ©. 

In the D 3> |/3| limit, Eq. ( |B5| ) peaks at A = with the value (3 2 /D 2 , with resonance 
width |AA| ~ D. This is to be compared with the familiar D = case in which the nominal 
resonance width |AA| ~ \/3\ corresponds to sin 2 2# m attaining a maximum value of unity at 
resonance. 
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FIGURES 



FIG. 1. Schematic representation of the polarisation P, the matter potential vector V, and the 
instantaneous diagonal basis S q in relation to the fixed coordinate system S. The authors warn 
that some features have been exaggerated for pictorial clarity. Note also that the unit vectors q\ 
and q\ are complex vectors. 

FIG. 2. Schematic representation of the precession axis (73 for the case \/3\ ~ D as a function 

of the matter potential vector V. The tags (/?, 0, A) show the pertinent values of V. 

FIG. 3. Schematic representation of the precession axis for the case D 3> \P\ as a function 

of the matter potential vector V. The tags (/?, 0, A) show the pertinent values of V. 

FIG. 4. Evolution of the effective total lepton number L", where a = fi,r, for v a <-> v s 
oscillation parameters Am 2 = — O.OleV 2 and sin 2 20$ = 10~ 7 . The dashed and dotted lines 
represent results from numerically integrating Eqs. ([34]) and (|35|) using the exact eigenvalue £3 
and the heuristically derived /c| tatlc respectively. These are juxtaposed with the solution to the 
exact QKEs [Eq. (Q)] for the same oscillation parameters (solid line). 

FIG. 5. Evolution of the effective total lepton number L^ a \ where a = fi, r, for u a *-* v s 
oscillation parameters Am 2 = — 10eV 2 and sin 2 2#o = 10~ 7 . The dashed and dotted lines represent 
results from numerically integrating Eqs. (^4|) and (|35| ) using the exact eigenvalue k% and the 
heuristically derived fc| tatlc respectively. These are juxtaposed with the solution to the exact QKEs 

[Eq. (0)] for the same oscillation parameters (solid line). 

FIG. 6. Damped harmonic oscillator with a "decaying" oscillation midpoint — schematic rep- 
resentation of the solution to Eq. (|59|). 

FIG. 7. Schematic representation of the matter and collision-affected mixing angle 2Q m ^. The 

ordinary matter- affected mixing angle is labelled 29 m . 

FIG. 8. The variable P y (t) as a function of A(i), where A(i) varies linearly with time from 

+A to —A, for adiabaticity parameters T = 0.25 (top left), 0.625 (top right), 1.25 (centre left), 

3.125 (centre right), 6.25 (bottom left), and 18.75 (bottom right). The parameters D and (3 are 

held constant at D = 0.02 and (3 = 0.001, and P z (t) — 1 throughout the evolution. The initial 

conditions are P x (0) ~ -Py(O) ~ 0, and P z (0) ~ 1. 
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